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I. INTRODUCTION 


The generalized circular error function, GCE, also known as the elliptical normal function, 
defines a probability function P(R,u,v). GCE gives the probability of a shot falling in a 
circle of the xy plane, of radius R and centered at the mean, under a bivariate uncorrelated 
normal distribution with mean zero and standard deviations u,v. This probability is given 
by 


PRuy =f eos pen GID? + 2) hava () 


Transforming to polar coordinates, with x = rcos@, y =rsin@ yields 


ff exp{—5 [SEE y? + (E82 rao ar. (2) 


Using elementary trigonometric identities 


PRY) 


27TUVv 


R/u 
P(R/u, c) -—| ie exp -5t °(B — Acos@)| rd 6dr, (3) 
with 298 bes 
_t-e Ae ree 
ar, Be: ” 


where, without loss of generality, v < u, 


c=v/u, Oc, (5) 
Using the fact that 
1 Tv 
Ip(x) = -{ exp(x cos 0) d 6, (6) 
T JO 


where Io(x) denotes the modified Bessel function of the first kind and zero order [1, p. 375], 
one obtains from (3) and (6) 


R/u 2 2 
PIR/) == f exp) (ran (7) 
C Jo 2 2 
Now letting 
K=R/u, war /2, (8) 
(7) reduces to 
1 K? /2 
PKG) = - | exp(—B w) Ip(A w) dw. (9) 
c Jo 


Introducing the MacLaurin series and the asymptotic expansion for Ig, [1, p.375, p.377] 
recurrence relations are developed which are used to generate an algorithm that efficiently 
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and accurately evaluates P (see [2], [8]). The objective of this report is to describe an 
algorithm upon which a Fortran computer program, INVGCE, is based to determine Kk, 
given P,u,v. This program can be used as a tool in evaluating the accuracy of a weapon. 

For example, consider a target located at the origin of the xy-plane and a set of miss 
distances normally distributed in the x and y directions (see [6], |7]) which are independent, 
with mean zero and standard deviations u, v respectively. Then, approximating the standard 
deviations from the data, say i and V, one is interested in knowing the radius, R, of the 
circle, C, centered at the origin, which contains 100 P% of the distribution, i.e., for which a 
shot has a probability P of falling within C. 

In Appendix A, the five-six significant digit inverse P(K,c) table in [2] is extended to 
six-seven significant digits. 


II. ANALYSIS TO DETERMINE K (R= Ku) 
INVGCE first treats two special cases. When c = 0, then P = erf(K//2), where 


erf(x) = a exp(—t?) dt, (10) 
T JO 
and if c = 1, then P = 1 — exp(—K?/2). 

For the first case, the inverse of the erf function is needed to determine K. The routine 
DERFI(erf(x), 1 — erf(x)) from [5, p.51] is used for this purpose. In the second case, K = 
\/—2log(1 — P). The remainder of this section assumes c € (0, 1). 

In order to find K in the general case, where c and P = P are given, the well-known 
Newton-Raphson root-finding procedure (N-R) is used [4, p. 129]. Denoting the nth approx- 
imation or iterate for K by K,, with c constant, 


K = =1,2,.-- 11 
n+1 Ka (OP /OK)x-x, ) n ioe) ) ( ) 
where if A K? < 14 (see [2]), then 
OP sai : 
OK — Cc exp( BK /2) I,(AK /2), (12) 
or else oP OK 
Fe =~ exp(-K?/2) [exp(—AK?/2) 1 (AK?/2)]. (13) 
c 


There are two basic problems associated with using N-R: 


1. A condition for stopping the iteration procedure is needed, i.e., deciding when an ac- 
ceptable approximation for K has been reached. 


2. A value for K,, which initiates the N-R, is needed so that rapid convergence is assured. 


NSWCDD/TR-04/43 


The first problem was resolved by requiring that the iteration be stopped when 


| Kati — Kn | < Ky+1 eps, (14) 


where eps is a prespecified small positive number. It was found unnecessary to utilize the 
stopping rule: 
| P(Kn1,c) — P| < Peps. (15) 


The second problem required extensive testing and some analysis to resolve satisfactorily. 
The solution adopted distinguishes this work from that carried out in [2] and [8]. The 
reasoning will be heuristic. 

Values for K, were considered using the results of [3], but those did not do as well, in 
general, as the values given below. 

From (9), the first term of the Taylor series for P is given by 

ZG 

142 


So, replacing T, by P, exp(—BK?/2) ~ T = 1-—- pe P, Ky, = /—2Log(T)/B. It was 
found, by extensive testing, that the expression for K, gave excellent starting values, provided 
T > .025. Thus 


‘Ts [1 — exp(BK?/2)]. (16) 


If T > .025, then K; = /(—2Log(T)/B). (IERR = 1). (17) 
If the inequality in (17) does not hold, the following set of conditions is used: 
If P < 10~° then 


K, = V2 DERFI(200 P, 1 — 200P). (IERR = 2) 
Else 
if P < 10-4 then 
K, = V2 DERFI(100P, 1 — 100P). (IERR. = 3) 
Else 
if P < 1073 then 
K, = V2 DERFI(60P, 1 — 60P). (IERR = 4) 
Else 
if P < 107! then 
K, = V2 DERFI(5P,1—5P). (IERR = 5) 
Else 
K, = V2 DERFI(P, 1 — P). (IERR = 6) 


IERR is an output parameter of INVGCE that identifies which K; was chosen. (See the 
next section.) The motivation for choosing K, in terms of the inverse of the erf function 
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arises from the fact that, for c = 0, the exact value for K is given in terms of this function, 
as noted by (10). 


III. DESCRIPTION OF THE FORTRAN 77 SUBROUTINE INVGCE 
In this section, the subroutine INVGCE is discussed further. Its call line is: 


Call INVGCE(P, c, K, ij, IERR) 


The routine, with 23 supporting routines, is written in Fortran 77 in double-precision. 
Thus, on an IBM PC approximately 15 decimal digits are available. The inputs are P and 
c. The outputs are K, ij, IERR. IERR plays a dual role as an error parameter and as an 
indicator of which initial approximation is used for Kj. 

If IERR is negative, then an unacceptable input has been given. K = —10!° is returned. 
In particular: 

If IERR=—-1 then P>1lorP <0. 
If IERR=-—-—2 then e< 0; 
If IERR=-—3 then CoA, 
The integer ij represents the number of calls to the N-R procedure. The output is K = R/u. 

INVGCE is designed to give a minimum of 8 significant digits for K when eps is set to 

5-10-°, and P is constrained to 
10°? P= 110. 


The routine is also efficient, showing, with extensive testing, a maximum for ij of 7 and an 
average of ij = 3. This indicates that the choices for K, are quite good. 
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APPENDIX A 


TABULATION OF K AS A FUNCTION OF P AND c. 
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TABLE A-1. TABULATION OF K AS A FUNCTION OF P and c 


Fixed value of P on a row. Fixed value of c on a column. 


Example 1: P=.31, c=.42, K=0.5800701 See Page A-13 
Example 2: P= .9990, c=.65, K = 3.3759637 See Page A-26 


A-3 
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